#include <complex>
#include <fstream>
#include <iostream>
#include <string.h>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<time.h>

using namespace std;


#define pi 3.1415926535898

template<typename T>inline T SQR(T x)
{  return(x*x);  }
template<typename T>inline T POW3(T x)
{  return(x*x*x);  }


template<typename T>inline int anint(T x)
{  return(int(floor(x+0.5)));  }

#include "random.h"
#include "procedimientos.h"

int main(void)
{
	/*Contadores*/
	int st, strel;
	/*numero total de particulas*/
	int npart;
	/*Posicion de la pared superior*/
	double ho,wd,def;
	/*velocidad de elongacion*/
	double vwd,vel;
	/*Incremento temporal*/
	double tfinal,trelax,deltat;
	/*polidispersidad*/
	double sigma,fi;
	double vol,rmax;
	/*posiciones de las particulas*/
	double **pos, *d;
	/*fuerzas*/
	double **force;

	npart = 100;
	ho=20.;
	tfinal=3000.;
	trelax=1500.;
	vel=0.01;
	sigma=0.0;
	fi=0.2;

	/*Creamos los arrays*/
	pos=new double*[npart];
	force=new double*[npart];
	d=new double[npart];
	for(int i=0;i<npart;i++){
		pos[i]=new double[3];
		force[i]=new double[3];
	}

	InitializeRandomNumberGenerator(1.0);
	/*Calculamos el diametro de las particulas*/
	/*(distribucion gaussiana media 1 variancia sigma)*/
	vol = setDiameters(npart, sigma, d);
	rmax=sqrt(vol/6./ho/fi);
	/*Colocamos las particulas*/
	setPositions(npart, pos, d, rmax, ho);
	/*Paso temporal*/
	deltat=5.0e-4;
	/*Numero de pasos temporales*/
	st=int(tfinal/deltat);
	/*Numero de pasos de relajacion antes de empezar la compresion*/
	strel=int(trelax/deltat);

	wd=ho;
	vwd=0.;
	/*Comienza a moverse la barrera*/
	for(int i=0;i<st;i++){
		if(i>strel){
			vwd=vel;
		} 
		/*Desplazamiento de la posicion de la barrera*/
		wd-=vwd*deltat;
		/*Calculo de las fuerzas entre particulas*/
		calc_force(npart, pos, d, force, wd, vwd);
		/*Desplazamiento de las particulas*/
		move_balls(npart, pos, force, d, deltat,wd);

		/*Almacenamiento de las posiciones de las particulas*/
		if(i%1000==0){
			time_t tiempo = time(0);
			struct tm *tlocal = localtime(&tiempo);
			char output[128];
			strftime(output,128,"%d/%m/%y %H:%M:%S",tlocal);
			def=(ho-wd)/ho;
			printf("%s: %5.5f	%5.5f\n",output,wd,def);
			calc_energy(npart, pos, d, ho, wd,i*deltat);
		}
	}

	/*Eliminamos los arrays*/
	for(int i=0;i<npart;i++){
		delete pos[i];
		delete force[i];
	}
	delete []pos;
	delete []force;
	delete []d;
}